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SUMMARY 

An efficient technique based on low-rank separated approximations is proposed for computation of three- 
dimensional integrals arising in the energy deposition model that describes ion-atomic collisions. Direct 
tensor-product quadrature requires grids of size 4000^ which is unacceptable. Moreover, several of such 
integrals have to be computed simultaneously for different values of parameters. To reduce the complexity, 
we use the structure of the integrand and apply numerical linear algebra techniques for the construction 
of low-rank approximation. The resulting algorithm is 10^ faster than spectral quadratures in spherical 
coordinates used in the original DEPOSIT code. The approach can be generalized to other multidimensional 
problems in physics. Copyright © 0000 John Wiley & Sons, Ltd. 
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1. INTRODUCTION 

Computation of multidimensional integrals is one of the most time-consuming tasks in physics. 
Standard approaches either have high complexity or require sophisticated quadrature schemes. 
Already in a three-dimensional case the integrand may depend on many parameters and should 
be computed many times, so the computational cost of the simplest tensor-product quadratures is an 
important issue. 

One of the promising tools to reduce the dimensionality of the problem (and hence the number 
of mesh points where the integrand should be calculated) is the usage of separation of variables 
in the integrand (see, for example, lH] |2)- The idea was known for a long time (for examples we 
refer the reader to the review |[3l), but it has become practically useful after the fast algorithms of 
decompositions of functions in a separated form had been obtained in two- aiiia, three- Eel and 
multidimensional cases EM- 

Let F{x,y) be a function of two variables x,y where point {x,y) is in a certain rectangle 
[ax,bx] C) [oy, by]. The function is said to be in a separated form if it can be represented as a sum of 
products of univariate functions; 

9 

F{T,y) =^arUr{x)gr{y). ( 1 ) 

r —1 

The minimal number q such that exists is called separation rank. Direct generalization 
of O to multivariate functions is referred to as a canonical polyadic (CP, also known as 
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CANDECOMP/PARAFAC) decomposition l|3][T0l. The reader can find examples of applications 
of separated representations in multidimensional cases in lfn fmfT^[T3l[T^fr5lfTHfT71fT8l[T9l . 

For a function given in the separated form the integration is simplified a lot. Indeed, 


P C pby 

/ / F{x, y)dxdy = E- Ur{x)dx / gr{y)dy, 
J J T—\ 


( 2 ) 


and the problem is reduced to the computation of one-dimensional integrals, which can be computed 
using fewer quadrature points than the original integral. 

In case of a discrete approximation of ([T]i 

9 

F{xi,yj) ^'^arUr{xi)gr{yj) ( 3 ) 

T = 1 


with the error e in the Frobenious norm the number q is called e-rank. We shall assume the notion 
of e-rank when the term rank will be mentioned in the text bellow. Expression (O can be rewritten 
in the matrix form 

A^UEG^, (4) 

where A is an n x m matrix with elements Ay = F{xi,yj), U is an n x q matrix with elements 
Uir = Ur{xi), G is an m X q matrix with elements Gjr = griyj) and T, is a q x q diagonal matrix 
with elements CTt- on the diagonal. This is a standard low-rank approximation problem for a 
given matrix. Provided that a good low-rank approximation exists, there are very efficient cross 
approximation algorithms aa that need only 0{(n + m)q) elements of a matrix to be computed. 

In this paper we describe how to apply this technique to speedup computations of three- 
dimensional integrals in the energy deposition model intended to describe ion-atomic collisions. 
This model was introduced by N. Bohr |[^ and developed further by A. Russek and J. Meli 1211, 
C.F. Cocke l22l, and V.P. Shevelko at al. Theoretical development of the model is presented 
in Il23ll24ll25ll^ . Examples of calculations are reported in Il27ll28ll29l[30l . Detailed description of 
the computer code DEPOSIT and user guide are given in ISTl and its updated version based on the 
separated representations is avialable in 1^ . 

The code DEPOSIT allows to calculate total and multiple electron loss cross sections and 
ionization probabilities needed for estimation of losses and lifetimes of fast ion beams, background 
pressures and pumping requirements in accelerators and storage rings. All of them are, in fact, 
functionals of the deposited energy T{h) (b is the impact parameter of the projectile ion), which in 
turn is a three-dimensional integral over the coordinate space. To calculate any of these parameters 
one has to compute T{b) in all points of the 6-mesh. 

In the original work ED an advanced quadrature technique was used, and the computational time 
has appeared to be much faster in a comparison with direct usage of uniform meshes. Calculation 
of the deposited energy T{b) for a given atomic shell in one point b took about several seconds. 
For complex ions it was about few hours on one processor core in total, that is not enough fast. 
To overcome this issue a fully scalable parallel variant of the algorithm was proposed, but the 
computational time was still large. 

We present an entirely different approach for the computation of T{b) in many points of the b- 
mesh based on the idea of separation of variables ([D. An approximation of functions to be integrated 
by a sum of products of univariate functions allows to effectively decrease the dimensionality of the 
problem. This requires active usage of numerical and analytical tools. 

The problem setting is as follows. The deposited energy T{b) is a three-dimensional integral 


T{b) 



z 


b)p{x, y, z)dxdydz. 


( 5 ) 


It involves cylindrically symmetric function of two variables (the energy gain Ai? during an ion- 
atomic collision) which is a smooth finite function and spherically symmetric function of three 
variables (electron density in a Slater-type approximation) which decays exponentially. For details 
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and definitions we refer the reader to Appendix [A] Previous approach ED used tensor-product 
quadratures in spherical coordinates. We use very fine uniform meshes and low-rank approximation. 

In Section ini the Slater-type function of three variables is decomposed by the exponential sums 
approach ll^ l3^ . The integral is immediately reduced to a two-dimensional one of a simpler 
structure. In Section lT^ for a function of two variables we use the pseudo-skeleton decomposition of 
matrices E5l|36l[33 computed via a variant of the incomplete cross approximation algorithm 
Combining these two representations we obtain then an efficient algorithm with 0{n) scaling, in 
comparison with 0{n^) complexity for direct integration over a three-dimensional mesh. We show 
numerically that the function in question can be well-approximated by a separable function. Thus, 
the approximation can be computed in 0{n) time, where n is the number of grid points in one 
dimension. The computation of T{h) on the whole &-mesh takes less then one minute instead of 
several hours and total speedup of the program is about ^ 10^ times. Illustrative examples are given 
in Section l23] All the equations related to the physical model are written in atomic units. 


2. NUMERICAL PROCEDURE 


2.1. Exponential sums. 

In this section we present analytical expansion of the spherically symmetric electron density p-yir) in 
Cartesian coordinates as a sum of separable functions. We use this decomposition later in Section l2^ 
for analytical integration in one dimension. 

A three-dimensional electron density Py(r) is taken in a Slater-type approximation 

Py{r) = r = +y‘^ + z'^, ( 6 ) 

where integer parameter and real Cy and fly correspond to one atomic shell labeled by 7 (see 
also Appendix |B] for description of parameters). Eor the density py{r) the following normalization 
condition occurs 

Py{r)dr = Ny, (7) 

where Ny is the number of electrons in a 7 -shell. Eor the sake of simplicity index 7 will be skipped 
and only one shell will be considered in equations bellow. 

Eor a function p{x, y, z) defined in (| 6 ll the separation of variables O in Cartesian coordinates can 
be done analytically Il3^l3^[^l39l . The main idea is to approximate the Slater density function by 
a sum of Gaussians 

K 

( 8 ) 

Once the approximation (O is computed, the separation of variables in Cartesian coordinates is 
immediately done 

K 

p{x, y,z)^Y^ Afc (9) 

fe =0 

The technique for computation of nodes A^ and weights pk is based on the computation of inverse 
Laplace transform. Let us consider a function fapit) such that its Laplace transform is a function 

Fapis): 

Fc^pis) = ^ e-^Va/3(i) dt, (10) 

where a and /3 are parameters of Slater density (| 6 ll. Inverse Laplace transform fap{x) can be 
computed analytically for known Fap{s) (fTOl) . In Appendix |B] we present explicit expressions for 
functions fapit) corresponding to function (fTOl i for integer a and real positive /3. 
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Once function fap {t) in expression (fTOl i is known, the integral (fTOl i can be computed numerically 
and approximated by a quadrature formula 

K 

p{r) Ki (11) 

k=0 

Here Wk and e**" are quadrature weights and nodes, respectively. The procedure to compute weights 
and nodes was taken from the paper |[34ll . For the reader’s convenience we give the formula and its 
derivation in Appendix ICl 

Comparision with equation (l8]l gives 

Afc = C'n;fce*Va/3(e*''), Vk = e**’. (12) 

It appears that not many quadrature points (at fixed r) are required to achieve the accuracy of the 
expansion of order 10“^ in the Chebyshev norm due to the hyper-exponential decay. Practically, 
such an accuracy is quite enough for the physical meaning of the model. However, r is going to 
be very small, as soon as the finer grids (i.e. for a large number of nodes) are required for higher 
precision. Estimation of upper bound K in sum (fTTI i which has the logarithmic dependence is given 
in Appendix Id] 


2.2. Fast computation of T{b). 

In this section we describe a numerical scheme based on the cross decomposition of two dimensional 
integrand. Dimensionality reduction (from three to two dimensions) is achieved by means of the 
separated representation of Slater density obtained in the previous Section. 

Discretization of one-dimensional integrals in (O by some quadrature formula with nodes Xi G 
[qx, bx], i = 1, ... yj G [oy, by], j = 1, ... ,m and weights leads to the approximation 


F{x, y)dxdy ; 


^ cr^ ^ ^ w'f\, 

r^l 2^1 j=l 


(%■)• 


(13) 


To get the decomposition Q we apply 2d-cross algorithm |l4]|40l implemented in ll32l . 

A three-dimensional integral T{b) defined in (|5]l can be reduced to a two-dimensional integral by 
means of the decomposition Q 


K 


T{b) = Y,^k 

k=0 



z 




(14) 


and analytical evaluation of one-dimensional Gaussian integral 



(15) 


K 


T{b) = v^Y. -^JJ ^E{x, z-h) dxdz. 


fc =0 


Suppose that AE{x, z — b) has been decomposed as follows 

<? 

AE{x, z — b) K. arUr{x)gr{z — b). 


(16) 


(17) 


T — 1 


Then the integration (fT^ can be reduced to a sequence of one-dimensional integrations. 


T{b) = V 4^ V arIrkJrkib), 


(18) 
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Irk-J 

2 

[ Urix)e dx, 

(19) 

Jrk (^) — J 

f griz — dz. 

du 

( 20 ) 

For the numerical approximation of integrals (fT9]l and (| 20 ]| we use quadrature formula with uniform 
quadrature nodes (although any suitable quadrature can be used) 

Irk = ^ W^^'''’ur(Xi)e 

i G O^(fc) 

= {f 1 > e} 

( 21 ) 

Xi — dx “t“ f' hx'^ 

Q<i<2Nr, hr=ax/Nx, 

( 22 ) 

Jrkib) = ^ w'f'’ QriZj - b)e , 

3 

(23) 

Zj = -a^+ j hz, 

0 < j < 2Nz, hz = az/Nz- 

(24) 

We sample the impact parameter b (which can take only positive values) with the same 

step hz 

bi=lh 

This allows us to introduce a new variable i 

z, 0<l<Nz. 

' = z — b discretized as follows 

(25) 

Zk = —2az + k hz, 0 < fc < 3Nz, 

and such that for the boundary conditions (l 2 ^. (l25Tl. (|26]| 

(26) 

- 

+ 

T 

jtT 

II 

"O 

(27) 


Thus, the approximation problem (fTTI i reduces to the low-rank approximation of the extended 

{2Nx + 1) X (3-/Vz + 1) matrix 


^ ^ (7-7-1/-,-(X; )g-7-(■^J ) ■ (28) 

T —1 

This should be done only once (using the cross approximation algorithm), and the final 
approximation of integral (|2^ reads 

Jrkibi) = ^ (29) 

j G nf(/c) 

The calculation of T{b) can be summarized in the following algorithm. 

1: for every 7-shell of the projectile ion do 
2: compute the decomposition (l8]l for p{r) 

3: compute the cross approximation for the matrix AE{xi,Zj) defined in (|28] | 

4 : for k = 0 .. .K do 

5: for T = 1 ... q do 

6: compute the integral Irk defined in 

7: for every bi required do 

8: for fc = 0 ... iT do 

9: for r = 1... g do 

10: compute the integral Jrkibi) defined in (l29l) 

11: compute T^{bi), equation (fTSl) 
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Finally, the question is how to compute T{b) for many different values of b. The direct summation 
requires 0{N^) operations. But a closer look reveals that it is in fact a discrete convolution, which 
can be computed in linear cost. Nevertheless, due to the exponential decay of factors and 
in sums and (|29] | correspondingly, there is a very small number of terms, such that a direct 
summation after the truncation is much faster. This can be easily seen from values of parameter 9^ 
defined as 


(27V, + l)(iT + l)’ 


(30) 


where N{X) is the cardinality of set X. For details, please, see Table |II] and discussion in the 
following Section l23] 


2.3. Numerical experiments and discussions. 

Once the analytical expansion (|9]l is obtained, the full calculation of T{b) consists of two steps (see 
algorithm in Section lT^ : calculation of the cross decomposition of integrand (fT^ on the extended 
z-mesh (|26] | and computation of all integrals Irk and Jrkibi) for all bi from (l25T l. There is no need in 
all the time recalculation of the cross approximation (l28T l for the full computation of these integrals. 
After it is computed (that is CTt-, Ur{xi) and grizj) are known), the calculation of Irk and Jrk{bi) 
starts. 

In Table H] we present times Tcross and ranks q calculated for the energy gains AEy{x,z) 
corresponding to two ion-atomic collision examples (Au^®+ + O and 17^®+ + Xe). Values of Tcross 
are reported for different sizes of {x, z)-mesh and the relative accuracy of the cross decomposition. 
One can find that Tcross is about 10“^ ~ 10° second in order of magnitude. Given that the 
full computation of integrals Irk and Jrk{bi) for all bi takes approximately 5000 • 5 • 10“® = 25 
seconds (the average value of column Ts in Table |II] for 5000 b values), we can conclude that 
the cross approximation becomes then a pre-computing stage with a tiny contribution to the total 
computational time. 

An important parameter in sum (|28] | is the rank value q. It determines the complexity of the 
algorithm (the smaller q, the better). As it follows from the numerical experiments, ranks of 
AEry{x, z) decomposition are small (Table Ul against the mesh size. It means that for real physical 
systems the cross decomposition is a very prominent tool. It allows to decrease the complexity of 
the problem in practice from 0{n^) elements to 0{q ■ n) elements where q<^n. 

In Table HI] we present the program speedup for every atomic shell. In quadrature sum 

all terms less then e = 10“^° were thrown out for every Xi due to the exponential decay of 
2 

factors for every fixed pk- Parameter 9^ is defined as a relative number of terms in total 

summation of Irk over all k and i (equation and d^l. above the threshold e. As it can be seen 
from the table, there are only a few percent of terms to be summed, which considerably reduces the 
sum and speeds up the total calculation. Parameter 9^ decreases when going from top to bottom of 
the third column (for one system). That is why the time Ts also decreases while both ranks K and 
q increase. The same situation occurs for sum d29] l over zj. 

Another important question is the full accuracy of the computation. As it was mentioned above, it 
consists of two contributions acquired from the cross approximation and the quadrature summations. 
In equations d2ni and d29l) values ar , Ur {xi ) and gr {zj ) are approximated with the cross accuracy Ec, 
while the quadrature sum is calculated with an error Si. In Table |III| and Table HVl we provide a 
numerical example demonstrating the actual error for the integral T{b) for small enough mesh sizes. 
We consider the worst case 6 = 0, when the integrand has the most singular behavior (because the 
Slater density has no higher derivatives at r = 0). As it follows from our results the scheme holds 
the third order for the Simpson rule, even in the worst singular case. In other cases 6 > 0 it holds the 
fourth order. These results show that the claimed accuracy is achieved and the numerical scheme 
holds the order of integration up to the accuracy of the approximated function. 

Finally, we can conclude that usage of the technique based on separated representations (fTSi) 
allows to decrease the total time of computing of T(6) by a factor of ^ 10® in comparison to the 
previous version. In practice, computational time is reduced from several hours to one minute or 
less on the same hardware. 
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Table I. Ranks q of the decomposition (I28l l calculated by incomplete cross approximation algorithm (4) 
for the energy gain AE{x,z). Two cases are considered: collision of ions with an Oxigen atom at 

a collision energy E = 6.5 MeV/u and collision of 17^®+ ions with a Xenon atom at a collision energy 
E = 2.5 MeV/u. The number of a;-mesh points are take n eq ual to 2N + 1, the number of z-mesh points 
is taken equal to ?>N + 1 in correspondence to equations (12211 and ( 126b . ax = az = 8. The relative error of 
the approximation in the Frobenius norm e and N are placed in the bottom line as a couple {£,N). They 
correspond to different numerical tests. Calculations were carried out on 1.3 GHz Intel Core i5 processor. 
Column Tcross corresponds to the time the cross algorithm takes. The numerical results confirm almost 

linear scaling of the approach in N. 


System 

7 -Shell 

q 

Tcross (sec) 

q 

Tcross (sec) 

q 

Tcross (sec) 

Au‘‘^+ + O 


13 

0.21 

21 

0.42 

24 

2.41 


4sp® 

13 

0.21 

21 

0.33 

24 

2.40 



14 

0.19 

22 

0.42 

26 

2.58 


3sp® 

16 

0.25 

24 

0.54 

29 

2.64 


2 sp® 

17 

0.28 

25 

0.56 

30 

2.71 


Isp'^ 

17 

0.27 

25 

0.56 

30 

2.70 

C/^®+ -h Ae 


14 

0.20 

22 

0.50 

26 

2.17 


4d/24 

15 

0.23 

24 

0.52 

27 

2.58 


4sp® 

17 

0.28 

25 

0.55 

30 

2.69 



17 

0.27 

25 

0.54 

30 

2.77 


3sp® 

17 

0.27 

25 

0.55 

30 

2.75 


2 sp® 

17 

0.26 

25 

0.54 

30 

2.71 


Isp'^ 

17 

0.26 

25 

0.55 

30 

2.69 



(10- 

1024) 

(10- 

-y, 1024) 

(10- 

4096) 


These results were obtained by using of our implementation of the cross approximation algorithm 
and the low-rank format of the Slater density (HUl- The latest version of the cross decomposition 
code implemented in C++ can be downloaded from the link 1401 . 


3. CONCLUSIONS AND FUTURE WORK 

A new and efficient technique for computation of three-dimensional integrals based on low-rank 
and separated representations is proposed for the energy deposition model. Due to a general form 
of the integrand, which is a product of two functions with cylindrical and spherical symmetries, 
this methodology can be applied to many types of integrals having similar structure. Such an 
approach significantly reduces computational time and allows to achieve a given accuracy (because 
it uses the cross approximation and quadrature summations). The general concept can be applied 
to more complicated models (like ion-molecular collisions with electron loss and charge-changing 
processes) that lead to multidimensional integrals. For the multidimensional case we plan to use fast 
approximation techniques based on the tensor train (TT) format ||9ll . 
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A. PHYSICAL MODEL 

In this section we introduce functions and parameters involved into the definition of deposited energy 
integral (O 

T^{h) = Ilh j{p)pj{x,y,z)dxdydz, (31) 
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Table II. Timings to compute T{b) at fixed b are presented for two cases: the DEPOSIT code (old) Td 
and the code based on separated representations ( fl8b Tg. Collision systems are the same as in Table U 
In the third column K labeles maximum value of index in expansion (O. The mesh ( 154b for radial 
density is taken with at = —3, bt — 45, ht = {bt — at)/250. The meaning of parameter Ox is explained 
in Section [Z3] Calculations were carried out for the relative accuracy of the cross decomposition e = 10“^ 
and [—8, 8] ® [—16, 8] mesh with 4097 x 6145 points. The last column shows speedup of the program. 


System 

7 -Shell 

K 

6x X 10 -^ 

Ts (xlO ® sec) 

Td (sec) 

TdITs 

+ O 


74 

9.1 

7.94 

3.89 

490 


4sp® 

69 

6.0 

4.92 

3.83 

778 



73 

4.3 

3.59 

3.88 

1080 



72 

3.9 

3.81 

3.82 

1003 


2sp® 

107 

1.5 

2.42 

3.86 

1592 


\sp^ 

209 

0.4 

1.24 

3.88 

3120 

+ Xe 

bsp* 

62 

13.1 

10.1 

3.94 

390 


4d/24 

70 

6.5 

6.05 

3.90 

644 


4sp® 

67 

5.1 

5.00 

3.94 

788 



71 

3.6 

3.88 

3.92 

1011 


isp^ 

70 

3.4 

3.52 

3.90 

1106 


2sp® 

105 

1.3 

1.99 

3.87 

1945 


\sp^ 

207 

0.3 

1.04 

3.88 

3723 


Table III. Convergence of integrals T{b) for 4d/^^, 4sp® and shells of + O at 6.5 MeV/u on two- 

dimensional mesh [—10,10[ (g) [—20,10[ with {2N + 1) x (3A^ + 1) points for different N (first column), 

see equations (122b and ( 126b . Simpson weights u;f, are used in (EB and ( 129b . For radial density the 
mesh ( 154b is taken with parameters at = —3, bt = 45, ht = {bt — at)/650. Calculations were carried out for 
the relative accuracy of the cross decomposition Sc = 10“^^ and fixed value of parameter feg = 0.0. Last 
column shows the rellative error . Extrapolated value of integral is calculated by Romberg method (with 
Richardson extrapolation) 1411 . p. 161. The order of scheme pe is defined by Aitken rule II41I . p. 344. 


7 -Shell 

N 

q 

Tibo) 

Pe 

Si 


2048 

30 

177.131769802401 


2.1 • lO-"^ 


4096 

33 

177.131804304375 


1.2- 10"® 


8192 

35 

177.131806364504 

4.07 

7.7- 10-1° 


16384 

37 

177.131806491972 

4.01 

4.8- 10-11 


32768 

39 

177.131806499918 

4.00 

3.0-10-12 


65536 

40 

177.131806500407 

4.02 

2.3- 10-1® 


extrapolated 

177.131806500448 



4sp® 

2048 

31 

165.507815905465 


2.2-10-1 


4096 

33 

165.507850781567 


1.3-10-® 


8192 

35 

165.507852865842 

4.06 

8.3- 10-1° 


16384 

37 

165.507852994821 

4.01 

5.2- 10-11 


32768 

39 

165.507853002863 

4.00 

3.2- 10-12 


65536 

40 

165.507853003364 

4.01 

2.1 - 10-1® 


extrapolated 

165.507853003399 




2048 

33 

407.589200892012 


5.0-10-1 


4096 

35 

407.589382823491 


5.2-10-® 


8192 

38 

407.589401536350 

3.28 

5.8-10-° 


16384 

40 

407.589403612947 

3.17 

6.7- 10-1° 


32768 

42 

407.589403853285 

3.11 

8.0-10-11 


65536 

44 

407.589403882005 

3.06 

9.8- 10-12 


131072 

46 

407.589403885498 

3.04 

1.2- 10-12 


262144 

47 

407.589403885925 

3.03 

1.9- 10-1® 


extrapolated 

407.589403886002 
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Table IV. Convergence of integrals T{b) for 3sp®, 2sp® and Is^ shells of + O at 6.5 MeV/u on two- 

dimensional mesh [—10,10] (g) [—20,10[ with {2N + 1) x (3V -I- 1) points for different N (first column), 

see equations (I22ll and ( 126b . Simpson weights wf, are used in (EB and ( 129b . For radial density the 
mesh ( 154b is taken with parameters at = —3, bt = 45, ht = {bt — at)/650. Calculations were carried out for 
the relative accuracy of the cross decomposition Sc = 10“^^ and hxed value of parameter feg = 0.0. Last 
column shows the rellative error . Extrapolated value of integral is calculated by Romberg method (with 
Richardson extrapolation) II41I . p. 161. The order of scheme pe is defined by Aitken rule 11411 . p. 344. 


7 -Shell 

N 

q 

nbo) 

Pe 

£i 

3sp^ 

2048 

35 

53.0906497215656 


3.0-10-'’ 


4096 

38 

53.0920135580097 


4.6-10"® 


8192 

41 

53.0922247714446 

2.69 

6.1 • 10 -'^ 


16384 

44 

53.0922529097596 

2.91 

7.7-10-® 


32768 

47 

53.0922564862273 

2.98 

9.7-10-® 


65536 

50 

53.0922569351057 

2.99 

1.2 • 10 -® 


131072 

52 

53.0922569912665 

3.00 

1.5- lO-''® 


262144 

54 

53.0922569982873 

3.00 

1.9- lO-''! 


524288 

56 

53.0922569991611 

3.01 

2.4- lO-''^ 


1048576 

58 

53.0922569992692 

3.01 

4.0- lO-''® 


extrapolated 

53.0922569992903 



2sp^ 

2048 

36 

4.97891259072213 


2.3-10-4 


4096 

39 

4.97986325621640 


3.5-10-® 


8192 

42 

4.98001545713255 

2.64 

4.7-10-® 


16384 

45 

4.98003568931593 

2.91 

5.9 -10-'^ 


32768 

48 

4.98003825660591 

2.98 

7.4- 10-® 


65536 

51 

4.98003857868963 

2.99 

9.2- 10-® 


131072 

53 

4.98003861898554 

3.00 

1 . 2 - 10 -® 


262144 

56 

4.98003862402352 

3.00 

1.4- 10-4° 


524288 

58 

4.98003862465325 

3.00 

1 . 8 - 10-44 


1048576 

60 

4.98003862473183 

3.00 

2.3- 10-42 


2097152 

61 

4.98003862474152 

3.02 

3.5-10-43 


extrapolated 

4.98003862474327 



Isp^ 

2048 

36 

0.122305706797573 


8.3 -10-3 


4096 

39 

0.123229108221677 


8.3-10-4 


8192 

42 

0.123322068141309 

3.31 

7.9-10-3 


16384 

45 

0.123330850494328 

3.40 

8 . 2 - 10 -® 


32768 

48 

0.123331746174083 

3.29 

9.1 - lO-^" 


65536 

51 

0.123331845477160 

3.17 

1.1 - 10 -'^ 


131072 

53 

0.123331857117543 

3.09 

1.3- 10-8 


262144 

56 

0.123331858525395 

3.05 

1 . 6 - 10 -® 


524288 

58 

0.123331858698474 

3.02 

2 . 0 - 10-4° 


1048576 

59 

0.123331858719929 

3.01 

2.5- 10-44 


2097152 

61 

0.123331858722600 

3.01 

3.1 - 10-42 


extrapolated 

0.123331858722980 




p = a/(2-6)2 _|_3,2_ (32) 

where the integration is done over the whole coordinate space. Integral 131b is written for an atomic shell 
with principal quantum number n and orbital quantum number I labeled by 7 = nh 

The energy gain AE-y{p) is a kinetic energy deposited to the projectile’s 7 -shell by the field U(R) of the 
target atom. 

AE.y(p) = AEj (p)nf(vj — 19) + AE^(p) nf(i9 — Vj), (33) 


nf(x) 


1 

g-kx _|_ I’ 


(34) 
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where fc = 3 by default and is an input parameter of the model. Expression ( 133b consists of two terms 
corresponding to low 

t 3 

AE<(p) = (35) 

P + C7^ 

and high 

<t>iFiXtP)^ , (36) 

energies. Function 

pOO - 

F^x) = xKi{x} = / dy (37) 

7o 

is defined in terms of the modified Bessel function of the second kind Ki{x) (see 142b . p. 375, Eq. 9.6.2) 
with asympthotical behavior (' ifel . p. 378, Eq. 9.8.3 and p. 379, Eq. 9.8.7) 

F{x 0+) = 1 + (i In I + 0.03860786) x^ + 0{x^). (38) 

Fitting parameters fi, Xi of the atomic field U{R) are obtained from the Dirac-Hartree-Fock-Slater 
calculations 1431 . The distance between the center of the field and the projectile electron of 7 -shell is labeled 
by R. Atomic field is given in the Yukawa potential form 

7 3 

UiR) = --J2<Pie~^''^, (39) 

2 = 1 


AS>(p) = 


P^ + P-y 


E 


where Z is the nuclear charge of the target atom. 

Parameter 1 ? is a relative velocity of the projectile. It is related with the collision energy of the ion-atomic 
system. The rest parameters v-^, u-y and /i-y concern to fundamental properties of the projectile ion 

and the target atom (such as binding energy, average velocity, mean radius of the shell, atomic radius and 
charge). They should be considered as positive constants in the integral (131b for a given ion-atomic system 
for all b. Detailed description of how to calculate them can be found in ED Examples of input files with the 
original code can be downloaded from link 1441 . 


B. INVERSE LAPLACE TRANSFORM SOURCES 


For integer a and real positive j3 inverse Laplace transform fapit) of F^p{s) from equation blOb may be 
calculated analytically and expressed via the Kummer’s confluent hypergeometric function M{a, b\ z) ( 1421 . 
chapter 13) as follows 


/a/3(f) 




(40) 


where 


M{a,b;z) = 1 -|- —Yy + 


a{a + 1 ) z^ 
b(b+l) ¥ 


-b ... 


(41) 


and r(a;) is the Gamma function. 

Below we present fap{t) explicitly for most practically usefull cases (a = 0,1,..., 6 ). Due to the 
difference of normalization conditions in spherical and Cartesian coordinates for the Slater density (lU 


p(r) = 


m^2rL-2pr 

r( 2 p + l) 


parameter a in ([^ is related to parameter p, in (142b as follows 


(42) 


a = 2p — 2. 


(43) 


The number of electrons in the shell 7 is labeled as N.y. Parameter p is greater or equal to unity. It is an 
integer or half-integer depending on the principal quantum number n and the orbital quantum number I 
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of the atomic shell. Details can be found in (Ulli. For example, = 1, a = 0; A‘ 2 sp 8 = 2, a = 2; 
/r 4 (^io = 3.5, O' = 5. Finally, 


For practical reasons higher 


2v^/33 ’ ^3/2 t) 

r u, 3^2 (V/3^) e-T 2\ 

^ 20F^4 ’ 92{t) ^5/2 3J 

4v^^5 ’ t5/2 i + 3f2j 

^ 15(?4(f//3^) e-7 4 , 4 N, 

40F/36 ’ 3t^l5t^) 

15g5{t/P^) e-i A 6 4 E 

80F^7 ’ S'sW ^7/2 (^1 ^+^2 45^3 

^ 105 56 (V/3') ,,, e-T 2 , 4 

’ 80F/38 ’ ^ i9/2 V t 5f2 105f3 

values of a are not necessarily due to the limitation of shell 


r) 


(46) 

(47) 

(48) 

(49) 


r) (50) 

filling with electrons. 


C. QUADRATURE FORMULA FOR THE LAPLACE INTEGRAL 
To obtain the decomposition ([ 8 ]l for given a and /3 we make a substitution s ^ into the equation JlOb 

f°° 2 

= / e“^ fo,j3{x)dx, (51) 

Jq 

then introduce another variable a: = e* in the integral 

/ CXD 2 t 

e”* " +*/a/ 3 (e*)df. (52) 

-OO 

The function under the integral ( 152b has exponential decay both in the spatial and frequency domains. 
Therefore the truncated trapezoidal rule gives the optimal convergence rate. It turns out to the final 
approximation of the form 

K 

Fap{s^) ~ ^ Wke’^fap{e'")e~‘‘ ® (53) 

k=0 


with trapezoidal weights Wk and the integrand values in the nodes e ® ® *^^**“/a/ 3 (e*''). The Gaussian-like 
part is split out in a separate factor in correspondence with decomposition ([ 8 ]l. Parameters of the formula 

tk = at + kht, ht = {bt-at)/K (54) 

have to be selected in such a way that the resulting quadrature formula approximates the integral for a wide 
range of parameter s. Typically, the choice at > —3, bt < 45, and K ~ 250 gives good accuracy (< 10“^). 
As an example, in Table |n] the required number of terms in sum ( 153b is presented. Accurate error analysis 
can be found in (Mi- 
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D. CONVERGENCE RATE OE K{N) EOR THE SEATER DENSITY SERIES 

In this section we estimate an upper bound K in sum (EB which has the logarithmic dependence on the 
mesh size N. We follow the proof given in paper im for Lemma 4.3. In this Lemma a Slater-type function 
p(-ys) from eq. dlOb is considered for a = 0. Below, this function is considered for integer a. 

The integrand in ( llOl l after the change of variables t = 

= (55) 

has the following decay on the real axis (skipping the numerical factor before the exponent) 

D ^ ^ -se^-cix , a-f 1 - (a mod 2) 

Pap{x)^e ", as a; 00 , ci =-^-, (56) 

Pcjp(a;) « as®—>-— 00 , C 2 = oi+^. (57) 

This immediately implies expression (5.3) from 1471 for b = min(/3^, s), h is taken in the notation of 1471 . 
Following then the statement I of Lemma 4.3 we may conclude, that K — 0{\ loge|(| log£| -f log 1/fe)) with 
the error e of the approximation and b ~ 1/A^. 
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